Référence pour les séries chronologique :
Hyndman, R.J., & Athanasopoulos, G. (2018) Forecasting: principles and practice, 2nd edition, OTexts: Melbourne, Australia. OTexts.com/fpp2. Accessed on 2 may 2022
Source : enquête fréquentation touristiques (https://www.ispf.pf/bases/Tourisme/EFT.aspx)
library(data.table)
library(ggplot2)
library(hts)
library(plotly)
data = fread("SerieTourisme.csv")
# la structure
str(data)
## Classes 'data.table' and 'data.frame': 8725 obs. of 4 variables:
## $ date : IDate, format: "2006-01-01" "2006-02-01" ...
## $ Region : chr "" "" "" "" ...
## $ PaysDiff : chr "" "" "" "" ...
## $ Touristes: int 108 56 47 46 54 24 21 41 21 34 ...
## - attr(*, ".internal.selfref")=<externalptr>
# quelques statistiques
data[,Region:=as.factor(fifelse(Region=="",NA_character_,Region))]
data[,PaysDiff:=as.factor(fifelse(PaysDiff=="",NA_character_,PaysDiff))]
summary(data)
## date Region PaysDiff
## Min. :2006-01-01 Europe (hors France):2993 Allemagne : 192
## 1st Qu.:2009-10-01 Asie :2113 Autre Europe: 192
## Median :2013-09-01 Pacifique :1336 Belgique : 192
## Mean :2013-09-26 Amérique du Sud : 742 Canada : 192
## 3rd Qu.:2017-08-01 Amérique du Nord : 573 France : 192
## Max. :2022-03-01 (Other) : 956 (Other) :7753
## NA's : 12 NA's : 12
## Touristes
## Min. : 1.0
## 1st Qu.: 14.0
## Median : 43.0
## Mean : 326.1
## 3rd Qu.: 194.0
## Max. :8748.0
##
# la série sous format TS
Touristes <- data[,lapply(.SD,sum,na.rm=T),.SDcols="Touristes", by="date"]
Touristes.ts <- ts(Touristes$Touristes,
frequency = 12,
start = Touristes[date==min(date), c(year(date),month(date))],
end = Touristes[date==max(date), c(year(date),month(date))])
print(Touristes.ts)
## Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2006 13714 15176 16256 18061 18051 18138 22917 20035 21144 21120 18401 18536
## 2007 15354 15859 18107 17290 17818 18783 21034 20897 18789 20891 16229 17190
## 2008 13102 14765 16829 15962 17485 16551 19111 18601 17989 18121 14040 13940
## 2009 9999 10372 12415 11230 13236 13824 16853 15808 14888 15972 12892 12958
## 2010 9016 9730 10547 10271 11525 12119 17790 15087 15160 16092 12784 13798
## 2011 11371 11038 12304 12458 12838 14424 16858 15372 14402 14519 13086 14106
## 2012 10238 11523 13075 13147 13879 14940 16979 16002 15944 15519 12470 15262
## 2013 11174 11177 13897 12011 13534 15120 17289 14655 14175 14576 12953 13832
## 2014 12422 12410 15410 15737 14853 14650 17656 14603 15500 17546 14646 15169
## 2015 12343 12949 14472 13956 14832 16223 18060 16463 16927 17561 15681 14364
## 2016 12340 14444 15168 17042 15738 16758 19517 17719 17117 17499 14295 14858
## 2017 11910 13564 16281 15523 16782 17599 21448 18563 18539 18271 15449 15030
## 2018 11457 15747 17452 15955 16559 19372 24168 20110 19809 20661 17241 17737
## 2019 15007 16752 18684 19240 18749 21487 25361 21864 20305 21170 19185 18838
## 2020 13948 15497 7491 4605 7834 7680 8976 4486 6500 3924 524 347
## 2021 3368 8552 14331 9481 6739 12346 10320 12321 6262 8982 14856 293
## 2022 13714 15176 16256
derAnnee <- floor(max(time(Touristes.ts)))
derMois <- round((max(time(Touristes.ts)) %% 1)*12 + 1)
ggplotly(
autoplot(Touristes.ts, ylab="", xlab="", main="Evolution du nombre de touristes") +
geom_point(data=window(Touristes.ts,c(derAnnee,derMois),c(derAnnee,derMois)), size=2,color="red",show.legend=F))
ggplotly(
ggseasonplot(Touristes.ts, ylab="", xlab="") +
ggtitle("Evolution du nombre de touristes par mois") +
guides(colour=guide_legend(title="Années",ncol=2)))
ggplotly(
ggmonthplot(Touristes.ts, ylab="", xlab="") +
ggtitle("Moyenne des variations annuelles par mois du nombre de touristes"))
FirstAnnee <- floor(min(time(Touristes.ts)))
FirstMois <- round((min(time(Touristes.ts)) %% 1)*12 + 1)
LastAnnee <- 2019
LastMois <- 12
Touristes.ts2 <- window(Touristes.ts,c(FirstAnnee,FirstMois),c(LastAnnee,LastMois))
ggplotly(
ggseasonplot(Touristes.ts2, ylab="", xlab="") +
ggtitle("Evolution du nombre de touristes par mois") +
guides(colour=guide_legend(title="Années",ncol=2)))
ggplotly(
ggmonthplot(Touristes.ts2, ylab="", xlab="") +
ggtitle("Moyenne des variations annuelles par mois du nombre de touristes"))
L’hypothèse retenue sur la saisonnalité est qu’elle ne varie pas au cours du temps. Le résidus est le résultat de la série d’origine à laquelle on soustrait la saisonnalité et la tendance estimées.
library(dplyr)
ggplotly(
Touristes.ts2 %>% decompose(type="additive") %>%
autoplot() + xlab("Année") +
ggtitle("Décomposition saisonnière classique")
)
Les fonctions à retenir :
* stl()
* seasonal()
* trendcycle()
* remainder()
* seasadj ()
Touristes.stl <- stl(Touristes.ts2, s.window="periodic", robust=TRUE)
remain.ts <- remainder(Touristes.stl)
remain.temps <- time(remain.ts)[which(abs(remain.ts)>quantile(abs(remain.ts),probs=0.95))]
remain.valeur <- remain.ts[which(abs(remain.ts)>quantile(abs(remain.ts),probs=0.95))]
remain.dt <- data.table(datetime=remain.temps,y=remain.valeur,parts=factor("remainder"))
ggplotly(
autoplot(Touristes.stl) + ggtitle("Décomposition de la série (data = seasonal + trend + remainder)") +
geom_point(data=remain.dt,aes(x=datetime, y=y), size=2,color="red",show.legend=F) +
xlab("") + ylab("Touristes")
)
La fonction seasadj peut prendre différents objets TS, comme la décomposition classique
ggplotly(
autoplot(cbind(brute=Touristes.ts2,CVS=seasadj(Touristes.stl)), ylab="Touristes") +
ggtitle("Série corrigée des variations saisonnières") +
scale_colour_manual(values = c("dark gray","blue")) +
xlab("")+
theme_bw() +
theme(legend.position = "bottom")
)
# les séries sous format TS
data[,Region2:=Region]
data[Region %in% c("France","Europe (hors France)"),Region2:="Europe"]
PaysImp <- data[,.(Touristes=sum(Touristes)),.(PaysDiff,Region2)][Touristes>=50000, .(PaysDiff,Region2)]
data[, c("Region3","PaysDiff2"):=.("Autre","Autre")]
data[PaysImp, c("Region3","PaysDiff2"):=.(i.Region2,i.PaysDiff), on=.(Region2,PaysDiff)]
data[,c("Region3","PaysDiff2"):=.(sprintf("%10s",abbreviate(Region3,10)),
sprintf("%7s",abbreviate(PaysDiff2,7)))]
print(data[,.(Touristes=sum(Touristes)),.(Region3,PaysDiff2)])
## Region3 PaysDiff2 Touristes
## 1: Autre Autre 370260
## 2: AmériqdNrd Canada 103476
## 3: AmériqdNrd USA 918050
## 4: Asie Japon 201141
## 5: Europe Allemgn 60058
## 6: Europe Espagne 51547
## 7: Europe Italie 142705
## 8: Europe Roym-Un 61146
## 9: Europe France 640486
## 10: Pacifique Austral 131777
## 11: Pacifique Nvll-Cl 62953
## 12: Pacifique Nvll-Zl 101660
TouristesParPays <- dcast(data, date ~ Region3 + PaysDiff2, fun.aggregate = sum, value.var = "Touristes",sep="")
TouristesParPays.ts <- ts(TouristesParPays[,-1],
frequency = 12,
start = TouristesParPays[date==min(date), c(year(date),month(date))],
end = TouristesParPays[date==max(date), c(year(date),month(date))])
ggplotly(autoplot(TouristesParPays.ts, ylab="", xlab="",colour = T,facets = T,
main="Evolution du nombre de touristes par pays",
space="free_y") + theme(strip.text.y = element_blank()))
TouristesParPays.hts <- hts(TouristesParPays.ts, characters = c(10,7))
str(TouristesParPays.hts)
## List of 3
## $ bts : Time-Series [1:195, 1:12] from 2006 to 2022: 1599 1724 1758 1528 1745 1630 1998 1796 2117 2030 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : NULL
## .. ..$ : chr [1:12] " Asie Japon" " Autre Autre" " Europe France" " Europe Italie" ...
## $ nodes :List of 2
## ..$ Level 1: int 5
## ..$ Level 2: 'table' int [1:5(1d)] 1 1 5 3 2
## .. ..- attr(*, "dimnames")=List of 1
## .. .. ..$ : chr [1:5] " Asie" " Autre" " Europe" " Pacifique" ...
## $ labels:List of 3
## ..$ Level 0: chr "Total"
## ..$ Level 1: chr [1:5] " Asie" " Autre" " Europe" " Pacifique" ...
## ..$ Level 2: chr [1:12] " Asie Japon" " Autre Autre" " Europe France" " Europe Italie" ...
## - attr(*, "class")= chr [1:2] "hts" "gts"
ggplotly(autoplot(aggts(TouristesParPays.hts,levels = 0), ylab="", xlab=""))
ggplotly(autoplot(aggts(TouristesParPays.hts,levels = 1), ylab="", xlab="",colour = T,facets = T) + theme(legend.position="none"))
TouristesParPays1.stl.list <- lapply(aggts(TouristesParPays.hts,levels = 1), FUN = function(x) stl(x, s.window = "periodic", robust = T))
TouristesParPays1.trend.list <- lapply(TouristesParPays1.stl.list, FUN = function(x) trendcycle(x))
ggplotly(autoplot(
ts(as.data.table(TouristesParPays1.trend.list),
frequency = 12,
start = TouristesParPays[date==min(date), c(year(date),month(date))],
end = TouristesParPays[date==max(date), c(year(date),month(date))]),
main="Tendance du nombre de touristes par région", ylab="", xlab=""))
Projections sur les 12 mois à venir issues d’un modèle ARIMA. Le modèle a été élaboré suite à l’analyse des graphes d’autocorrélation et d’autocorrélation partielle. Le modèle a été testé sur plusieurs échantillons tests et mis en compétition avec un modèle de sélection automatique des paramètres.
# Autocorrélation et Autocorrélation partiel
Touristes.ts2 %>% ggtsdisplay(main="")
# Autocorrélation et Autocorrélation partiel sur série désaisonnalisée
Touristes.ts2 %>% diff(lag=12) %>% ggtsdisplay(main="Série désaisonnalisée")
# Autocorrélation et Autocorrélation partiel sur série désaisonnalisée et sans tendance
Touristes.ts2 %>% diff(lag=12) %>% diff() %>% ggtsdisplay(main="Série désaisonnalisée et sans tendance")
# Modèle arima
fitarima <- arima(Touristes.ts2, order=c(0,1,4), seasonal = list(order=c(0,1,1), period=12))
checkresiduals(fitarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,4)(0,1,1)[12]
## Q* = 19.803, df = 19, p-value = 0.4065
##
## Model df: 5. Total lags used: 24
# AIC pour comparer les modèles
AIC(fitarima)
## [1] 2586.577
# Projection
fc <- forecast(fitarima, h=12 , level=c(95))
autoplot(fc) +
ggtitle("Projection du nombre de touristes selon un modèle ARIMA",
subtitle = "Paramètres du modèle p=0, d=1, q=4, P=0, D=1, Q=1, s=12") +
guides(fill=guide_legend(title = "Intervalle \nde confiance" )) +
ylab("") +
xlab("")
Projections sur les 12 mois à venir issues d’un modèle hierarchique à 2 niveaux (pays et régions). Pour chaque série de la hierarchie les projections de base sont obtenues à l’aide d’un modèle automatique ARIMA. Les projections de base sont révisées de telle sorte que les aggrégations par niveau soient conséquentes. Par exemple, l’aggrégat des projections des pays d’une région correspond aux projections de la région en question.
library(parallel)
TouristesParPays2.hts <- hts(window(TouristesParPays.ts,c(FirstAnnee,FirstMois),c(LastAnnee,LastMois)) , characters = c(10,7))
fc.hts <- forecast(TouristesParPays2.hts, h=12, fmethod="arima", parallel = T, num.cores=detectCores())
ggplotly(
autoplot(allts(fc.hts, forecasts = T), ylab="", xlab="")
)